library(tidyverse)
library(phyloseq)
library(speedyseq)
library(ggrepel)
library(ampvis2)
library(plotly)
library(microbiome)
library(here)
options(getClass.msg=FALSE) # https://github.com/epurdom/clusterExperiment/issues/66
#this fixes an error message that pops up because the class 'Annotated' is defined in two different packages
rm(list = ls())
'%!in%' <- function(x,y)!('%in%'(x,y))
source("https://raw.githubusercontent.com/fconstancias/DivComAnalyses/master/R/phyloseq_taxa_tests.R")
source("https://raw.githubusercontent.com/fconstancias/DivComAnalyses/master/R/phyloseq_normalisation.R")
source("https://raw.githubusercontent.com/fconstancias/DivComAnalyses/master/R/phyloseq_alpha.R")
source("https://raw.githubusercontent.com/fconstancias/DivComAnalyses/master/R/phyloseq_beta.R")
source("https://raw.githubusercontent.com/fconstancias/DivComAnalyses/master/R/phyloseq_heatmap.R")
ps = "data/raw/metabarcoding/ps_silva_dada2_human_chicken_meta.RDS"
ps %>%
here::here() %>%
readRDS() %>%
phyloseq_get_strains_fast() %>%
filter_taxa(function(x) sum(x > 0) > 0, TRUE) %>%
subset_samples(Enrichment == "NotEnriched") %>%
subset_samples(Model == "Human") %>%
subset_samples(Treatment %!in% c("positive", "negative", "DONOR", "CCUG59168", "HV292.1", "STRAIN")) -> physeq
## Joining, by = "ASV"
physeq %>%
sample_data() %>%
data.frame() %>%
DT::datatable()
physeq %>%
rarefy_even_depth(sample.size = 4576,
rngseed = 123) -> physeq_rare
## `set.seed(123)` was used to initialize repeatable random subsampling.
## Please record this for your records so others can reproduce.
## Try `set.seed(123); .Random.seed` for the full vector
## ...
## 38 samples removedbecause they contained fewer reads than `sample.size`.
## Up to first five removed samples are:
## IR-13-S103IR-46-S126IR-51-S157IR-74-S188IR-82-S143
## ...
## 738OTUs were removed because they are no longer
## present in any sample after random subsampling
## ...
physeq_rare
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 372 taxa and 123 samples ]:
## sample_data() Sample Data: [ 123 samples by 59 sample variables ]:
## tax_table() Taxonomy Table: [ 372 taxa by 8 taxonomic ranks ]:
## phy_tree() Phylogenetic Tree: [ 372 tips and 371 internal nodes ]:
## refseq() DNAStringSet : [ 372 reference sequences ]
## taxa are rows
physeq_rare %>%
phyloseq_compute_bdiv(phylo = FALSE, norm = "pc") -> bdiv_list
## Loading required package: ape
## Loading required package: GUniFrac
## Loading required package: vegan
## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.5-7
##
## Attaching package: 'vegan'
## The following object is masked from 'package:microbiome':
##
## diversity
## Loading required package: matrixStats
##
## Attaching package: 'matrixStats'
## The following object is masked from 'package:dplyr':
##
## count
physeq_rare %>%
phyloseq_plot_bdiv(bdiv_list,
m = "CoDa", # with this potion I am generating CLR transformed OTU table
seed = 123) -> coda
bdiv_list$coda <- coda$physeq_clr %>% phyloseq::distance(method = "euclidean") # Aitchinson distance wich claim to be superior to the other ones for compositional data.
coda$PCA$layers[[1]] = NULL
coda$PCA + geom_point(size=2,
aes(color =Treatment,
fill = Treatment,
shape = NULL,
alpha = Day_of_Treatment)) +
theme_light() +
geom_path(arrow = arrow(type = "open", angle = 30, length = unit(0.15, "inches")),
size = 0.05, linetype = "dashed", inherit.aes = TRUE, aes(group=Reactor_Treatment, color = Treatment), show.legend = FALSE) +
scale_alpha_continuous(range=c( 0.9, 0.3)) + scale_color_viridis_d(na.value = "black") +
scale_fill_viridis_d(na.value = "black") +
# scale_shape_manual(values = c(8, 21, 22, 23, 24, 16, 15, 18, 17)) +
theme_classic() +
labs(col=NULL, fill = NULL, shape = NULL) + guides(shape=FALSE) -> pca_atchi
pca_atchi
plotly::ggplotly(pca_atchi)
# pca_atchi %>%
# export::graph2ppt(append = TRUE,
# file = file.path(here::here("data/processed/figures_NRP72")))
This is the fonction of interest
phyloseq_add_taxa_vector(dist = bdiv_list$coda,
phyloseq = physeq_rare,
figure_ord = pca_atchi,
pval_cutoff = 0.05,
top_r = 16, # top features tp display(strain, or depending on what you specified based on correlation r value)
taxrank_glom = "Strain",
tax_rank_plot = "Strain"
) -> out
Add the vectors to the ordination plot
out$plot
Get the raw table of taxa and correlation with the ordination:
out$envfit %>%
DT::datatable()
You could then filter that table and use the ASV id to generate a heatmap of only the significant features:
out$envfit %>%
filter(pval <= 0.05) %>%
pull(id) -> ASV_envit
physeq_rare %>%
prune_taxa(x = ., taxa = ASV_envit) %>%
phyloseq_ampvis_heatmap(transform = FALSE,
group_by = "SampleID",
facet_by = c("Model", "Fermentation", "Reactor_Treatment"),
tax_aggregate = "Genus",
tax_add = NULL,
ntax = 50) -> p
p$data %>%
mutate(Abundance = na_if(Abundance, 0)) -> p$data
p + facet_grid( ~Model + Fermentation + Reactor_Treatment , scales = "free", space = "free") +
scale_fill_viridis_c(breaks = c(0, 0.01, 0.1, 0.2, 0.50, 0.75, 0.100),
labels = c(0, 0.01, 0.1, 0.2, 0.50, 0.75, 0.100),
trans = scales::pseudo_log_trans(sigma = 0.001),
na.value = 'transparent') -> p2
## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.
p2
# p2 %>%
# export::graph2ppt(append = TRUE,
# file = file.path(here::here("data/processed/figures_NRP72")))
vegan::adonis(bdiv_list$coda ~ get_variable(physeq_rare, "Reactor_Treatment") *
get_variable(physeq_rare, "Day_of_Treatment"))$aov.tab %>%
data.frame() %>%
rownames_to_column('Term') %>%
mutate_if(is.numeric, round, 4) %>%
DT::datatable()
lapply(
bdiv_list,
FUN = phyloseq_adonis_strata_perm,
physeq = physeq_rare,
formula = paste0(c("Reactor_Treatment", "Day_of_Treatment"), collapse=" * "),
nrep = 999,
strata = "none"
) %>%
bind_rows(.id = "Distance") %>%
mutate_if(is.numeric, round, 3) %>%
DT::datatable()
phyloseq_plot_bdiv(physeq_rare,
bdiv_list,
m = "NMDS",
axis1 = 1,
axis2 = 2) -> plots
# plot_list %>%
# phyloseq_ordinations_expl_var()
physeq_rare %>%
phyloseq_distance_boxplot(bdiv_list$coda ,
d = "Reactor_Treatment") -> dist_box
dist_box_plot <- dist_box$plot
dist_box_plot
#
# dist_box_plot %>%
# export::graph2ppt(append = TRUE,
# file = file.path(here::here("data/processed/figures_NRP72")))
plots %>%
phyloseq_plot_ordinations_facet(color_group = "Day_of_Treatment",
shape_group = "Treatment") -> plot_list
plot_list
lapply(
bdiv_list,
FUN = physeq_pairwise_permanovas,
physeq = physeq_rare,
compare_header = "Reactor_Treatment",
n_perm = 999,
strat = FALSE
) %>%
bind_rows(.id = "Distance") %>%
mutate_if(is.numeric, round, 3) %>%
# filter(! terms %in% (c("Residuals", "Total"))) %>%
DT::datatable()
correct_plot_x_continious <- function(original_plot = tmp_plot$data){
ggplot(data=original_plot,mapping=aes(x=as.numeric(as.character(varGroup2)),y=Distance,color=Label)) +
# geom_boxplot(data=tmp_plot$data,mapping=aes(x=as.factor(varGroup2),y=Distance,color=Label),outlier.size = 0.5) +
geom_point(position=position_jitterdodge(jitter.width = 0.1,seed=123),aes(group=Label),size=2, alpha=0.4) +
theme_bw() +
geom_path(arrow = arrow(type = "open", angle = 30, length = unit(0.1, "inches")),
size = 0.4, linetype = "dashed", inherit.aes = TRUE, aes(group=Label, color = Label),
position=position_jitterdodge(dodge.width=0.9)) +
theme(legend.position="bottom") + geom_smooth(show.legend = FALSE, level = 0.95, alpha=0.005, size = 0.001) +
# ylim(c(0,1)) +
xlab("Day Treatment") + guides(col = guide_legend(ncol = 3)) + scale_x_continuous(breaks=seq(0,90,10)) + scale_fill_viridis_d() + scale_color_viridis_d() -> plot
return(plot)
}
phyloseq_plot_beta_div_wrt_timepoint(distances = c("coda", "bjaccard", "wjaccard"),
bdiv_list = bdiv_list,
physeq = physeq_rare,
timepoint = "fixed",
group_var = "Reactor_Treatment",
time_var = "Day_of_Treatment",
fixed_time = -2) -> t
tmp_plot <- t$bjaccard + facet_null()
tmp_plot$data %>%
correct_plot_x_continious(.) -> tmp_plot_corrected
tmp_plot_corrected
# tmp_plot_corrected %>%
# export::graph2ppt(append = TRUE,
# file = file.path(here::here("data/processed/figures_NRP72")))
tmp_plot <- t$wjaccard + facet_null()
tmp_plot$data %>%
correct_plot_x_continious(.) -> tmp_plot_corrected
tmp_plot_corrected
#
# tmp_plot_corrected %>%
# export::graph2ppt(append = TRUE,
# file = file.path(here::here("data/processed/figures_NRP72")))
tmp_plot <- t$coda + facet_null()
tmp_plot$data %>%
correct_plot_x_continious(.) -> tmp_plot_corrected
tmp_plot_corrected %>%
plotly::ggplotly()
# tmp_plot_corrected %>%
# export::graph2ppt(append = TRUE,
# file = file.path(here::here("data/processed/figures_NRP72")))
phyloseq_plot_beta_div_wrt_timepoint(distances = c("coda","bjaccard", "wjaccard"),
bdiv_list = bdiv_list,
physeq = physeq_rare,
timepoint = "previous",
group_var = "Reactor_Treatment",
time_var = "Day_of_Treatment") -> t
tmp_plot <- t$bjaccard + facet_null()
tmp_plot$data %>%
correct_plot_x_continious(.) -> tmp_plot_corrected
tmp_plot_corrected
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : span too small. fewer data values than degrees of freedom.
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -1.04
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 5.04
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 9.2416
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : span too small. fewer
## data values than degrees of freedom.
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## -1.04
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius 5.04
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 0
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 9.2416
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 2
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 2
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : Chernobyl! trL>n 6
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : Chernobyl! trL>n 6
## Warning in sqrt(sum.squares/one.delta): NaNs produced
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at 2
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius 2
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 0
## Warning in stats::qt(level/2 + 0.5, pred$df): NaNs produced
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf
tmp_plot_corrected %>%
export::graph2ppt(append = TRUE,
file = file.path(here::here("data/processed/figures_NRP72")))
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : span too small. fewer data values than degrees of freedom.
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -1.04
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 5.04
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 9.2416
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : span too small. fewer
## data values than degrees of freedom.
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## -1.04
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius 5.04
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 0
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 9.2416
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 2
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 2
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : Chernobyl! trL>n 6
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : Chernobyl! trL>n 6
## Warning in sqrt(sum.squares/one.delta): NaNs produced
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at 2
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius 2
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 0
## Warning in stats::qt(level/2 + 0.5, pred$df): NaNs produced
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf
tmp_plot <- t$wjaccard + facet_null()
tmp_plot$data %>%
correct_plot_x_continious(.) -> tmp_plot_corrected
tmp_plot_corrected
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : span too small. fewer data values than degrees of freedom.
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -1.04
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 5.04
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 9.2416
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : span too small. fewer
## data values than degrees of freedom.
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## -1.04
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius 5.04
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 0
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 9.2416
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 2
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 2
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : Chernobyl! trL>n 6
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : Chernobyl! trL>n 6
## Warning in sqrt(sum.squares/one.delta): NaNs produced
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at 2
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius 2
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 0
## Warning in stats::qt(level/2 + 0.5, pred$df): NaNs produced
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf
tmp_plot_corrected %>%
export::graph2ppt(append = TRUE,
file = file.path(here::here("data/processed/figures_NRP72")))
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : span too small. fewer data values than degrees of freedom.
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -1.04
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 5.04
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 9.2416
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : span too small. fewer
## data values than degrees of freedom.
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## -1.04
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius 5.04
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 0
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 9.2416
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 2
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 2
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : Chernobyl! trL>n 6
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : Chernobyl! trL>n 6
## Warning in sqrt(sum.squares/one.delta): NaNs produced
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at 2
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius 2
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 0
## Warning in stats::qt(level/2 + 0.5, pred$df): NaNs produced
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf
tmp_plot <- t$coda + facet_null()
tmp_plot$data %>%
correct_plot_x_continious(.) -> tmp_plot_corrected
tmp_plot_corrected
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : span too small. fewer data values than degrees of freedom.
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -1.04
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 5.04
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 9.2416
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : span too small. fewer
## data values than degrees of freedom.
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## -1.04
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius 5.04
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 0
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 9.2416
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 2
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 2
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : Chernobyl! trL>n 6
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : Chernobyl! trL>n 6
## Warning in sqrt(sum.squares/one.delta): NaNs produced
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at 2
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius 2
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 0
## Warning in stats::qt(level/2 + 0.5, pred$df): NaNs produced
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf
tmp_plot_corrected %>%
export::graph2ppt(append = TRUE,
file = file.path(here::here("data/processed/figures_NRP72")))
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : span too small. fewer data values than degrees of freedom.
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -1.04
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 5.04
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 9.2416
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : span too small. fewer
## data values than degrees of freedom.
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## -1.04
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius 5.04
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 0
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 9.2416
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 2
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 2
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : Chernobyl! trL>n 6
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : Chernobyl! trL>n 6
## Warning in sqrt(sum.squares/one.delta): NaNs produced
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at 2
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius 2
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 0
## Warning in stats::qt(level/2 + 0.5, pred$df): NaNs produced
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf
phyloseq_plot_beta_div_wrt_timepoint(distances = c("coda","bjaccard", "wjaccard"),
bdiv_list = bdiv_list,
physeq = physeq_rare,
timepoint = "between.ref.group",
group_var = "Reactor_Treatment",
time_var = "Day_of_Treatment",
group_to_compare = "CR_UNTREATED") -> t
tmp_plot<- t$bjaccard + facet_null()
tmp_plot$data %>%
correct_plot_x_continious(.) + guides(col = guide_legend(ncol = 2)) -> tmp_plot_corrected
tmp_plot_corrected
tmp_plot_corrected %>%
export::graph2ppt(append = TRUE,
file = file.path(here::here("data/processed/figures_NRP72")))
tmp_plot <- t$wjaccard + facet_null()
tmp_plot$data %>%
correct_plot_x_continious(.) + guides(col = guide_legend(ncol = 2)) -> tmp_plot_corrected
tmp_plot_corrected
tmp_plot_corrected %>%
export::graph2ppt(append = TRUE,
file = file.path(here::here("data/processed/figures_NRP72")))
tmp_plot <- t$coda + facet_null()
tmp_plot$data %>%
correct_plot_x_continious(.) + guides(col = guide_legend(ncol = 2)) -> tmp_plot_corrected
tmp_plot_corrected
tmp_plot_corrected %>%
export::graph2ppt(append = TRUE,
file = file.path(here::here("data/processed/figures_NRP72")))
sessionInfo()
## R version 4.0.2 (2020-06-22)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS High Sierra 10.13.6
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] gdtools_0.2.2 usedist_0.4.0.9000 GUniFrac_1.1
## [4] matrixStats_0.58.0 vegan_2.5-7 lattice_0.20-41
## [7] permute_0.9-5 ape_5.4-1 reshape2_1.4.4
## [10] scales_1.1.1 here_1.0.1 microbiome_1.10.0
## [13] plotly_4.9.2.1 ampvis2_2.6.5 ggrepel_0.8.2
## [16] speedyseq_0.5.0 phyloseq_1.34.0 forcats_0.5.0
## [19] stringr_1.4.0 dplyr_1.0.4 purrr_0.3.4
## [22] readr_1.4.0 tidyr_1.1.2 tibble_3.0.6
## [25] ggplot2_3.3.3 tidyverse_1.3.0.9000
##
## loaded via a namespace (and not attached):
## [1] readxl_1.3.1 uuid_0.1-4 backports_1.2.1
## [4] systemfonts_0.3.2 plyr_1.8.6 igraph_1.2.6
## [7] lazyeval_0.2.2 splines_4.0.2 crosstalk_1.1.0.1
## [10] digest_0.6.27 foreach_1.5.1 htmltools_0.5.1.1
## [13] magrittr_2.0.1 cluster_2.1.0 doParallel_1.0.16
## [16] openxlsx_4.2.3 Biostrings_2.56.0 modelr_0.1.8
## [19] officer_0.3.14 prettyunits_1.1.1 colorspace_2.0-0
## [22] blob_1.2.1 rvest_0.3.6 haven_2.3.1
## [25] xfun_0.21 crayon_1.4.1 jsonlite_1.7.2
## [28] survival_3.2-3 iterators_1.0.13 glue_1.4.2
## [31] rvg_0.2.5 gtable_0.3.0 zlibbioc_1.34.0
## [34] XVector_0.28.0 webshot_0.5.2 Rhdf5lib_1.10.1
## [37] BiocGenerics_0.34.0 DBI_1.1.1 miniUI_0.1.1.1
## [40] Rcpp_1.0.6 xtable_1.8-4 viridisLite_0.3.0
## [43] progress_1.2.2 stats4_4.0.2 DT_0.15
## [46] truncnorm_1.0-8 htmlwidgets_1.5.3 httr_1.4.2
## [49] RColorBrewer_1.1-2 ellipsis_0.3.1 pkgconfig_2.0.3
## [52] NADA_1.6-1.1 farver_2.0.3 dbplyr_1.4.4
## [55] manipulateWidget_0.10.1 later_1.1.0.1 tidyselect_1.1.0
## [58] labeling_0.4.2 rlang_0.4.10 munsell_0.5.0
## [61] cellranger_1.1.0 tools_4.0.2 cli_2.3.0
## [64] generics_0.1.0 devEMF_4.0-2 ade4_1.7-16
## [67] export_0.3.0 broom_0.7.2 fastmap_1.1.0
## [70] evaluate_0.14 biomformat_1.7.0 yaml_2.2.1
## [73] knitr_1.31 fs_1.5.0 zip_2.1.1
## [76] rgl_0.100.54 nlme_3.1-149 mime_0.10
## [79] xml2_1.3.2 compiler_4.0.2 rstudioapi_0.13
## [82] zCompositions_1.3.4 reprex_0.3.0 stringi_1.5.3
## [85] highr_0.8 stargazer_5.2.2 Matrix_1.2-18
## [88] multtest_2.44.0 vctrs_0.3.6 pillar_1.4.7
## [91] lifecycle_1.0.0 data.table_1.13.6 flextable_0.5.11
## [94] httpuv_1.5.4 patchwork_1.1.0 R6_2.5.0
## [97] promises_1.1.1 network_1.16.0 IRanges_2.22.2
## [100] codetools_0.2-16 ggnet_0.1.0 MASS_7.3-52
## [103] assertthat_0.2.1 rhdf5_2.32.4 rprojroot_2.0.2
## [106] withr_2.4.1 S4Vectors_0.26.1 mgcv_1.8-32
## [109] parallel_4.0.2 hms_1.0.0 grid_4.0.2
## [112] rmarkdown_2.4 Cairo_1.5-12.2 Rtsne_0.15
## [115] shiny_1.5.0 Biobase_2.50.0 lubridate_1.7.9
## [118] base64enc_0.1-3